#/usr/bin/perl
use POSIX "fmod"; # float modulus
use List::Util qw[min]; # min function
# GLOBAL
my $grid_size_x = 300;
my $grid_size_y = 300;
my $locitot = 99;
my $name = "sim";
my $vertical_scaling_factor = 1.5; # for hex grid is 0.75, we'll take modulus of 1.5
my $biopsy_size_x = 20; # size of biopsy, 10x10 = 100 cells maximum
my $biopsy_size_y = 10; # biopsies are 1mm by 2mm correspond to 10 crypts by 20 crypts


use Math::Geometry::Planar;
sub round {
    my($number) = shift;
    return int($number + .5 * ($number <=> 0));
}

open(out,">cell.xy.coordinates.per.distance.class.csv");
# Draw without replacement 100 cells from 1, 2, .., 300 cells apart = total 30,000 pairs

for $distance (1..300){

    my $cell_count = 0;
    my %cellxy = ();
    while($cell_count<100){
	
	# Pick a cell
	$rx = int(rand(150));
	$ry = int(rand(300));

	# If cell was already picked
	if($cellxy{$rx}{$ry}){
	    next;
	}

	# Center of circle
	$center = [$rx,$ry];
	$radius = $distance;
	# Compute coords of 100 cells that lie on the cirlce
	$ncells = 100;
	$polygon = CircleToPoly($ncells,$center,$radius);
	@points = $polygon->points();

	# Pick 1 cell at random from circle
	$circle_cell = int(rand($ncells));
	$rx2 = round($points[0][$circle_cell][0]);
	$ry2 = round($points[0][$circle_cell][1]);

	# Check if cell is within bounds
	# If it is, success, put 2 cells in hash
	if($rx2>=0 && $rx2<150 && $ry2>=0 && $ry2<300 && !$cell{$rx2}{$ry2}){	    
	    $cell_count++;
	    $cellxy{$rx}{$ry} = 1;
	    $cellxy{$rx2}{$ry2} = 1;
	    # randomly shift x, because of the wrap-around issue, make it more even on x-axis
	    $xshift = int(rand($grid_size_x));
	    # rx = 149, xshift = 299,  448-299, wrap around
	    if($rx+$xshift<299){$rx+=$xshift;} else {$rx=$rx+$xshift-299;}
	    if($rx2+$xshift<299){$rx2+=$xshift;} else {$rx2=$rx2+$xshift-299;}
	    #$xshift = int(rand(min($grid_size_x-$rx, $grid_size_x-$rx2)));	    
	    my ($hx, $hy) = abs_to_hex_coords($rx, $ry);
	    my ($hx2, $hy2) = abs_to_hex_coords($rx2, $ry2);
	    print out "$distance,$hx,$hy,$hx2,$hy2\n";
	    	    
	}
    }
}

# convert coords
sub abs_to_hex_coords{
    my ($x, $y) = @_;    
    $y = $y*($vertical_scaling_factor/2.0);
    # Depending on picking y that is divisible by 1.5 or not, add 0.5 to x                                         
    if(fmod($y,$vertical_scaling_factor)){
	$x = $x + 0.5;
    }
    return ($x, $y);
}
